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Abstract The Two Point Exponential Approximation Method was introduced 
by Fadel et al. (Fadel, 1990), and tested on structural optimization problems 
with stress and displacement constraints. The results reported in earlier 
papers were promising, and the method, which consists in correcting Taylor 
series approximations using previous design history, is tested in the present 
paper on optimization problems with frequency constraints. The aim of the 
research is to verify the robustness and speed of convergence of the Two 
Point Exponential Approximation method when highly non-linear constraints 
are used. 

Introduction 

In the practice of optimization, especially when complex structural, thermal, 
aerodynamic or other analyses are needed, the computer time required to 
perform the analyses is critical. Most large optimization problems have been 
formulated such that the number of full scale analyses are minimal. This is 
generally accomplished by reducing the original problem to an approximate, 
simpler model which can be optimized within certain constraints. The 
original problem is then solved with the optimized approximate design 
variables, and iterations are performed until overall convergence is attained. 
The critical aspect of the procedure is the quality of the approximation. For a 
very highly non-linear problem, linear approximations are valid only in a 
very small domain around the original design point, whereas in better 
behaved problems, larger moves can be accomplished. The trade-off 
between the quality of approximation and number of real analyses is what 
dictates the overall time needed for reaching the optimum (if at all 
reachable). 

Derivation of the Two Point Exponential Approximation 

Several traditional approximation methods were summarized in the paper by 
Fadel et al (Fadel, 1990) ranging from the simple Taylor series in the form: 
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g(X) = g(Xo) + X (X i- Xoi) 

i dxi 


to the reciprocal, hybrid, and higher order approximations. The authors then 
introduced the Two Point Exponential approximation which is an extension of 
the simpler Taylor series, adjusted by matching the derivatives at the 
previous design point. This correction term is incorporated into an exponent 
which is computed after each real analysis for each constraint, and with 
respect to each design variable. The exponent acts as a measure of goodness 
of fit: If the linear approximation is valid for a certain constraint, the 

exponent is close to or equal to 1, if the reciprocal approximation is more 
appropriate, the exponent approaches or is equal to -1. In other cases, the 
exponent varies between -1 and 1, correcting the approximation and 
improving the fit of the data. 

The Two Point Exponential Approximation is derived as mentioned earlier by 
matching the slopes at previous design points. Initially, one substitutes x Pi 
for x in the Taylor series: 


g(X) = g(Xo) + £ (xj*- x£j) 

i 3xf 


and after resubstitution, one can write: 


i 


Xfii 

Pi 


dg(X 0 ) 

dxj 


with the exponent evaluated according to: 


Pi = 


( 

(ag(Xi)U 

log 

l a** i 

[ag(Xo)| 

1 

\ axi ]) 



The point Xi refers to the design point at the previous iteration and X 0 refers 
to the current design point from where the approximation is carried out. 

Note that at the first iteration, since no previous design history exists, a 
linear or reciprocal step is carried out, depending on the preference of the 
user. 
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The results reported in the earlier paper compared the linear, reciprocal and 
Two Point Exponential approximations on structural problems with stress 
and displacement constraints. Three problems of different sizes were used, 
namely the standard three bar truss problem, a 25 bar truss transmission 
tower, and a 52 bar truss tower. The results showed that the Two Point 
Exponential approximation generally displayed a much smoother behavior 
than the other two methods. It contributed to reducing the oscillations 
between successive iterations, and required less iterations to reach the 
optimum in most cases. The overhead involved in computing the exponents 
proved to be insignificant. The exponents have to be computed after each 
real analysis and used during the optimization of the approximate problem. 
Care has to be taken in the code development to avoid divisions by zero, and 
to avoid having to compute the logarithm of a negative number. In such 
cases, the algorithm should be written in a way that it would revert to a 
linear or reciprocal step. 

Two Point Exponential Approximation and Frequency Constraints 

After ascertaining the merit of the approximation in the case of stress and 
displacement constraints, it was suggested to test the method on frequency 
type constraints. The frequency constraints are generally highly non-linear, 
and further testing of the method was warranted to confirm its value for 
general structural optimization problems. For this purpose, two test 
problems of different complexity and size were selected. The approximation 
method is tested on both the problems, and results and conclusions are 
reported. Both problems were taken from the literature to ensure 
correctness. 

Test problem: Cantilever Beam 

The first test problem is taken from Pritchard and Adelman (Pritchard, 

1990). The 193 inch long hollow cantilever beam with square cross section 
(Figure 1) has four design variables: the height and width of the beam cross- 
section, and the two wall thicknesses (sides, top and bottom). The beam is 
divided into ten elements. The first element near the base has a slightly 
different modulus of elasticity, but all other characteristics are uniform over 
the length of the beam. The dimensions and physical characteristics of the 
standard beam Xq are: 


H = 5.00 in 

B = 3.75 in 

t = 0.80 in 

d = 0.10 in 
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and moduli of elasticity: (Element 1 is at the wall) 

E 2-10 = 5.85E6 

Ei = 4.90E6 



The problem was analyzed using the ANSYS (Swanson, 1990) finite element 
package, and optimizations were carried out with the program CONMIN 
(Vanderplaats, 1973). The first test consisted in evaluating the 
approximations to the first bending frequency when one variable was 
modified. The height of the beam cross section was selected as the design 
variable, and the results are tabulated in Table 1. 


H 

linear 

reciprocal 

exact 

2 pt exp. 

H 

rel err lin [%] 

rel err 2pt [%] 

2 

1 .7303 

-3.64094 

1.67448 

1.62987 

2 

1.050261 153 

0.83990199 

3 

2.9239 

1.332393 

2.88837 

2.88399 

3 

0.668228188 

0.08239831 

4 

4.1175 

3.81906 

4.10689 

4.10836 

4 

0.199018652 

0.027637706 

5 

5.3111 

5.31106 

5.31106 

5.31106 

5 

0 

0 

6 

6.5047 

6.305727 

6.49736 

6.49678 

6 

0.137449021 

0.010981808 

7 

7.6983 

7.016203 

7.66538 

7.66855 

7 

0.619085456 

0.059730265 

8 

8.8919 

7.54906 

8.81562 

8.82852 

8 

1.435494986 

0.242964664 

9 

10.085 

7.963504 

9.94882 

9.97827 

9 

2.572744424 

0.554566968 

1 0 

11.279 

8.29506 

1 1.0658 

11.119 

1 0 

4.01539429 

1.002004366 


dfdh0= 1 

(H=5) 

1.1936 



P= 

0.929378792 


dfdhl = 1 

(H-6) 

1.17833 






Table 1. Cantilever Beam analysis. Evaluation of approximations 

based on Beam Height H. 
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These results are illustrated in Figure 2. The errors resulting from both the 
linear and Two Point Exponential approximations are plotted as a function of 
the beam height H. The reference point Xo is the point H=5, and the Two 
point Exponential approximation uses the previous analysis point at Xi as 
the point where H=6. The graphs show the superior performance of the new 
approximation in the case of changes in one design variable. 



Figure 2. Approximation error as a function of beam height H 

Cantilever Beam Problem 


When one considers variations in multiple variables simultaneously, the 
advantages of one approximation versus another are less easily 
demonstrated. In this study, we first considered changing all four variables 
of the cantilever beam problem simultaneously by progressive percentages. 
Figure 3 illustrates the relative errors of the approximations of the first 
bending frequency as a function of the relative change in all four design 
variables. These changes are certainly not indicative of performance within 
an optimization exercise, but they do provide some measure of goodness 
easily displayable. Note that the starting point for the approximation is the 
point with abscissa 0. The linear approximation is carried out from this point 
forward (increasing all four design variables by x%), and then backward 
(decreasing all four variables by x%). For the Two Point Exponential 
approximation, the starting point is the same Xo, and the "previous" design 
point Xi is at abscissa x=.2. In this case, increasing x means backtracking, 
whereas decreasing x means progressing in the direction established by the 
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two successive design points. The figure shows that the Two Point 
Exponential approximation seems to better fit the real function below the 
design point, and is slightly worse than the linear approximation above the 
design point. It is hoped that during an optimization, the design variables 
would either increase or decrease monotonically, and the Two Point 
Exponential approximation would perform better than the linear. Note that 
the results of both approximations were very sensitive to the derivatives 
obtained through finite differences in the analysis program (ANSYS). 

The true test of an approximation however, is to perform the optimization 
exercise. This is the subject of the next section. 



Deviation in design points (H, B, T, D) 


Figure 3. Approximation error as a function of all four variables 

Cantilever Beam Problem 


Optimization of the Cantilever Beam Problem. 

Since a true test of the approximations can only be obtained in an 
optimization problem, the Cantilever Beam example discussed above was 
reformulated as an optimization exercise. The initial design variables are the 
ones given above as vector Xo, and the object of the problem is to find the 
minimal weight subject to frequency constraints. The first frequency 
constraint is the first bending frequency of the beam which has to be below a 
certain minimum value and the second frequency above another value. This 
would ensure a separation of natural frequencies, and could be used as a 
design problem. The first attempt to solve the problem considered two 
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design variables, namely the height and width of the beam, leaving the 
thicknesses constant. The constraints (first and second frequencies) are 
limited to 5 Hz and 30 Hz respectively (FI < 5Hz, F2 >30Hz). The allowable 
error is 0.1 and the move limits are 50% in all three cases. The results are 
tabulated below: 



Linear 

Reciprocal 

2 Pt. Exp. 

0 

6.68 

6.68 

6.68 

1 

3.60146 

3.62224 

3.60146 

2 

2.13858 

2.16792 

2.15689 

3 

1.51095 

1.56556 

1.56387 

4 

1 .83129 

1.55081 

1.5507 

5 

1.52809 

1.55081 

1.5507 

6 

7 

1 .54845 



Table 

2. Variation 

of Cross 

sectional area 

as function of 

iteration 

number. 


Because of the similarity of results, a graph of the variation of objective 
(cross sectional area) with respect to iteration number would not provide any 
additional information. From the table above, one can only deduce that in 
this particular case, the three approximations perform relatively similarly. 

All three reach the optimum in roughly the same amount of steps. The linear 
approximation seems to reach a smaller optimum, but this result is because 
this particular approximation in this problem causes one of the constraints to 
be slightly violated, and at the final result, the second frequency constraint is 
active, but very close to be violated, whereas in the two other methods, the 
second frequency constraint is active, and satisfied. Table 3 lists some of the 
results for the above problem. In all three cases, the beam width is driven to 
the minimum (0.5 in), and the second frequency constraint becomes active. 



Linear B 

Linear F2 

Reciprocal B 

Reciprocal F2 

2 Pt Exp B 

2 Pt Exp F2 

1 

3.75 

33.6109 

3.75 

33.6109 

3.75 

33.6109 

2 

1.875 

29.6286 

1.875 

30.3661 

1.875 

29.6286 

3 

0.9375 

29.1512 

0.9375 

30.0812 

0.9375 

29.7322 

4 

0.5 

28.9012 

0.5 

30.41 18 

0.5 

30.3652 

5 

0.75 

28.1872 

0.5 

30.005 

0.5 

30.0019 

6 

0.5 

29.3767 

0.5 

30.005 

0.5 

30.0019 

7 

0.5 

29.9399 






Table 3 Cantilever Beam. Active constraints as function of iteration number. 
Beam width B driven to >= 0.5 in, second natural frequency driven to >= 30Hz. 


When one considers all four parameters: height, width and thicknesses, as 
design variables, the problem should be more complicated and the 
approximations less well behaved. 
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** Linear Reciprocal ” 2 Pi. Exp. 

Figure 4. Cantilever beam. Cross sectional area (objective) 
Four design variables, two constraints. 


Figure 4. illustrates the variation of the objective with respect to iteration 
number in the case of four design variables. In this case again, all three 
approximations behave relatively similarly, with the linear and Two Point 
Exponential approximations reaching the minimum in 5 steps whereas the 
reciprocal this time, takes one additional step. The interesting observation 
however, is the path to the minimum taken by all three approximations. In 
order to show the differences, the value axis (area) was magnified with a 
maximum at 2 inches. The first two iterations are therefore not visible, but 
one can see that the Two Point Exponential approximation is the smoothest 
behaved function. 

Figure 5. illustrates in the same problem (four design variables, two 
constraints), the variation of the second natural frequency. The problem 
consisted in minimizing the area subject to the second frequency remaining 
above 30Hz. The figure shows that the Two Point Exponential method shows 
similar oscillative behavior as the other methods, but with a smaller 
amplitude. 

The two results described sofar show that for two relatively simple problems 
with frequency constraints, the Two Point Exponential approximation 
behaves at least as good, if not better than the best of the linear or reciprocal 
approximations. 
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F 



Linear Reciprocal “ 2 Pt. Exp. 


Figure 5. Cantilever beam problem, four design variables 
Variation of second frequency constraint 


Conclusion 

The Two point Exponential Approximation was tested on problems with 
frequency constraints. The results obtained sofar show that the method is at 
least as performing as the best of the traditional methods like the linear or 
reciprocal approximation. It does also perform as a more controlled method 
which should be used when the problem to be solved does not have 
uniformly linearly behaved or uniformly reciprocally behaved constraints 
and objectives. 
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APPENDIX A 



Numerical 

Results 

of Optimization 

runs 

Linear approximation 
H B 

OBJ 

<-5 

CONI 

>=30 

COM2 

0 

5 

3.75 

6.68 

5.31106 

33.6109 

1 

4.60731 

1.875 

3.60146 

4.68022 

29.6286 

2 

4.79292 

0.9375 

2.13858 

4.60463 

29.1512 

3 

5.15476 

0.5 

1.51095 

4.56506 

28.9012 

4 

4.75644 

0.75 

1.83129 

4.45203 

28.1872 

5 

5.24046 

0.5 

1.52809 

4.64033 

29.3767 

6 

7 

5.34226 

0.5 

1.54845 

4.7295 

29.9399 

Reciprocal Approximation 
0 5 3.75 

6.68 

5.31106 

33.6109 

1 

4.71122 

1.875 

3.62224 

4.797 

30.3661 

2 

4.93959 

0.9375 

2.16792 

4.75187 

30.0812 

3 

5.4278 

0.5 

1.56556 

4.80423 

30.41 18 

4 

5.35406 

0.5 

1.55081 

4.73982 

30.005 

5 

6 

5.35406 

0.5 

1.55081 

4.73982 

30.005 

Two Point Exponential Approximation 
0 5 3.75 6.68 

5.31106 

33.6109 

1 

4.60731 

1.875 

3.60146 

4.68022 

29.6286 

2 

4.88447 

0.9375 

2.15689 

4.69663 

29.7322 

3 

5.41935 

0.5 

1.56387 

4.79686 

30.3652 

4 

5.35349 

0.5 

1.5507 

4.73932 

30.0019 

5 

6 
7 

5.35349 

0.5 

1.5507 

4.73932 

30.0019 


Cantilever Beam results 2 variables 2 constraints 









Linear 


<=5 

>-30 


H 

B 

T 

D 

OBJ 

CONI 

COM2 

1 

5 

3.75 

0.8 

0.1 

6.68 

5.31106 

33.6109 

2 

4.25085 

1.875 

0.4 

0.05 

1.84509 

4.69586 

29.7274 

3 

4.3655 

0.9375 

0.2 

0.05 

0.77155 

4.37265 

27.6857 

4 

5.27703 

0.5 

0.1 

0.05 

0.6077 

4.46111 

28.2445 

5 

5.65837 

0.5 

0.1 

0.05 

0.64584 

4.75696 

30.1133 

6 

7 

5.65837 

0.5 

0.1 

0.05 

0.64584 

4.75696 

30.1133 

1 

5 

3.75 

0.8 

Reciprocal 

0.1 

6.68 

5.31106 

33.6109 

2 

4.13298 

1.875 

0.4 

0.05 

1.8333 

4.56184 

28.8809 

3 

4.41934 

0.9375 

0.2 

0.05 

0.77693 

4.42169 

27.9955 

4 

6.34962 

0.5 

0.1 

0.05 

0.71496 

5.2913 

33.4862 

5 

5.70696 

0.5 

0.1 

0.05 

0.6507 

4.7946 

30.351 

6 

5.63605 

0.5 

0.1 

0.05 

0.64361 

4.73967 

30.0041 

7 

8 

5.63605 

0.5 

0.1 

0.05 

0.64361 

4.73967 

30.0041 

1 

5 

3.75 

0.8 

2 Pt Exponential 

0.1 6.68 

5.31 106 

33.6109 

2 

4.25085 

1.875 

0.4 

0.05 

1.84509 

4.69586 

29.7274 

3 

4.45606 

0.9375 

0.2 

0.05 

0.78061 

4.45508 

28.2064 

4 

5.84612 

0.521375 

0.1 

0.05 

0.66889 

4.92356 

31.1652 

5 

5.63885 

0.5 

0.1 

0.05 

0.64389 

4.74184 

30.0178 

6 

5.63885 

0.5 

0.1 

0.05 

0.64389 

4.74184 

30.0178 


Cantilever Beam results. 4 variables 2 constraints 
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APPENDIX B 

Ansys input file and Program listing 


H=5 . 

B=3 .75 
T=0 . 8 
D=0 . 1 

/TITLE, Beam model for approximation testing 3D model 

FINISH 

/PREP7 

KAN, 2 

KAY, 1,-1 

KAY, 2, 3 

KAY, 7, 3 

C*** compute area and IZZ 
IYY1=(T**3) *B 
IYY2=IYY1/12 
PAR1=H- (T*2 ) 

PAR3=B-(D*2) 

PAR 2= (H-T ) 12 
IYY3=(T*B) *(PAR2**2) 

IYY4=(IYY2+IYY3) *2 
IYY5=D*(PAR1**3) 

IYY6=IYY5/6 
IYY =IYY6+IYY4 
AREA=( (T*B)+(PAR1*D) ) *2 
C*** end of calculations 

ET,1,3 * 2D elastic beam 

R, 1, AREA, IYY, H 

HP, EX, 1,4. 9e6 * material properties for element 1 

MP,DENS, 1,0. 00018 

MP,EX,2,5.85e6 * material properties for other elements 

MP,DENS, 2, 0.00018 

N, 1, 0 

N, 11, 193 

FILL 

/PNUM, NODE , 1 

NPLOT 

MAT, 1 

E, 1, 2 

MAT, 2 

E, 2 , 3 

EGEN ,9,1,2 
EPLOT 
D, 1, ALL 

M,2,UY, 11,UX,R0TZ 
SAVE 

ITER, 1,1 

SFWRITE 

FINISH 

/SOLVE 

FINISH 
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/P0ST1 
set, , 1 

♦get, frel, freq 
set , , 2 

♦get, fre2, freq 
set, , 3 

♦get, fre3, freq 

FINISH 

/OPT 

FACT=. 99999 

H1=H*FACT 

H2=H/FACT 

B1=B*FACT 

B2=B/FACT 

OPVAR, H,DV, HI, H2 

OPVAR , B , DV , B1 , B2 

OPVAR, AREA, OBJ 

OPVAR, PARI , SV, . 1 , H 

OPVAR, PAR3 , SV, . 1 , B 

OPVAR, FRE1,SV, .1,10 

OPVAR, FRE2 , SV, . 1 , 100 

OPVAR, FRE2 , SV, . 1 , 150 

OPCOPY 

H=H*1. 001 

RUN, 2 

B=B*1 . 001 

H=H/ 1.001 

RUN, 3 

T=T*1. 001 

B=B/1. 001 

RUN, 4 

D=D*1. 001 

T=T/ 1.001 

RUN, 5 

OPLI ST , ALL , , 1 

FINISH 

/EOF 



c 

C23456789012 3 4 5678 9012 3 4 56789 012 3456789012345678901234567890123456789012 
C123456 7 

C 

C program to read an ANSYS file and extract the necessary data for 

C optimization, call conmin, and use approx to solve approximate prob 

c 

C Georges Fadel Sept 1990 

C Oct 1990 

C Jan 1991 

C 


C 


C 

C 


C 


C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
commons for CONMIN call 

COMMON /CNMN1/ DELFUN , DABFUN , FDCH , FDCHM , CT , CTMIN , CTL , CTLMIN , 

1 ALPHAX , ABOB J 1 , THETA , OBJ , NDV , NCON , NSIDE , IPRINT , NFDG , NSCAL , LINOB J , 

2 ITMAX, ITRM, ICNDIR, IGOTO,NAC, INFO , INFOG , ITER 

COMMON /CNMN2/ X (6) , DF (6) ,G (15) , ISC(15) , IC(15) , A(6, 15) , AF (7) 

COMMON /CNMN4/ VLB(6) ,VUB(6) ,SCAL(6) 

the next two are for approx subroutine. Second common just to pass 
flags to approx 

COMMON /INFOIN/ DV(4) , FUNC(7) , GRAD(4,7) 

COMMON /INFOLD/ DV1(4), FUNC1(7), GRAD1(4,7) 

COMMON /FLAGS/ IFLAG, ICALL, IDEBUG 
DIMENSION P (4 , 7) , RATX(4) , RATDER(4,7) 

DIMENSION S (6) ,G1(15) ,G2(15) ,B(15,15) ,C(15) ,MS1(30) 
end of conmin non-executable 
DIMENSION CONS (5, 5) ,OOBJ(4) ,GMAX(6) 

CHARACTER* 4 START (5) ,T (5) 

CHARACTER* 1 2 FILNM, FILNM1 , FILNM2 , FILNM3 
CHARACTER* 80 TT 
LOGICAL TOF 

DATA START ( 1) , START (2) , START (3) , START (4) , START (5) / 'LIST' , ' OPT', 

1 'IMIZ' , 'ATIO' , 'N SE'/ 

name of file (File=' ') written from batch file into 

temp.dat is read into FNAME. 

some parameters that have to be set for each optimization program: 
nlines in output file 
number of design variables NDV 
Number of constraints NCON 

Increment factor used to compute finite differences in 
finite element program: FACT = 1. - actual FACT 
DF means derivative of objective wrt design variable 
A means derivative of constraint wrt design variable 
and remember to adjust dimensions to read all needed data 
in X (NDV) , CONS (NCON, NDV) , OOBJ(NDV) 

DF (NDV) , A (NDV, NCON) 

Also, the output data includes a maximum of 6 cases per row. If 
NCON is more than 6, then, an additional read statement has to t 
written for the next batch of results. 

conmin requirements +++++++++++++++++++++++++++++ 

IGOTO Sets start of optimization loop 
IPRINT Print control: 0 print nothing 

1 print initial and final function informat 

2 1st debug level print 1 + control paramet 
function value and X at each iteration. 

3 2nd debug level print 2 + constraints, ac 
or violated constraints, move parameters. 



c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


approaches 0 as optimum gets closer 
4 full debug 

NDV Number of decision variables 

ITMAX Max number of iterations 

NCON Number of constraint functions G(J) 

NSIDE Number of side constraints (upper, lower bounds) 

G Constraints at initial design point 

CONS Constraints at finite differences 
ICNDIR Conjugate direction restart parameter 
NSCAL Scaling control parameter 

NFDG Gradient calculation control parameter 0: calculated by F 

Is externally supp 
2: obj external, r 

FDCH Relative change of decision variable for FD calc. 

FDCHM Minimum step for FD 
CT Constraint thickness parameter 

CTMIN Minimum abs value of CT 

CTL Constraint thickness for linear and side constraints 

CTLMIN Minimum abs value of CTL 

THETA Mean value of push off factor (for highly non-linear probl 
NACMX1 Estimate of number of active constraints 
DELFUN Minimum change in OBJ to indicate convergence 
DABFUN Same as DELFUN, but absolute not relative error 
LINOBJ 0 means non-linear, 1 means linear 

ITRM (3) number of consecutive iterations for convergence 

X(N1) Vector of decision variables 
VLB (Nl) Lower bound on variables X(I) 

VUB(Nl) Upper bound on variables X(I) 

SCAL(N5) Vector of scaling parameters not used if NSCAL=0 
ISC(N2) Linear constraint identification vector 

GMAX(NCON) LIMITS OF CONSTRAINTS 

IPRINT=2 SUPPLIED IN EXTERNAL FILE 
NDV=4 SUPPLIED IN EXTERNAL FILE 

SET NUMBER OF CONSTRAINTS TO REQUIRED NUMBER (INITIALLY 1, THEN 6) 
NC0N=1 SUPPLIED IN EXTERNAL FILE 


IGOTO=0 
NFDG=0 
ITMAX=50 
C NACMX1=15 
NSIDE=8 
ICNDIR=0 
NSCAL=0 
LINOBJ=0 
Nl=6 
N2=15 
N3=15 
N4=15 
N5=30 
ITRM=3 
FDCH=0 . 
FDCHM=0 . 

CT=0 . 

CTMIN=0 . 
CTL=0. 
CTLMIN=0 . 
THETA=0. 
DELFUN=10 . E-8 



DABFUN=10.E-8 

NAC=0 

ALPHAX=0 . 1 
ABOBJ1=0 . 1 
ICALL=1 
DO 9 1=1, N2 
ISC(I)=0 
9 CONTINUE 

C end of conmin variables definition 
C ++++++++++++++++++ 

NLINES=20000 
FACT=0 .001 

C ++++++++++++++++++ MAKE SURE THIS IS CORRECT TOTAL NUMBER OF 

C AND INCREMENT FACTOR IN ANSYS FOR DERIVATIV 

C create a file called OPTIM.DAT in which the filenames of the 

C initial result data file and file to be used to store results are 

C written. One name on each line. Next, enter a number representing 

C the magnitude of the move limits in % 

C 

OPEN (UNIT=3 , STATUS='OLD' , F I LE=' 0PTIM.DAT' ) 
read (3,99) FILNM , FILNM1 , FILNM2 , FILNM3 
r ead ( 3 , 9 8 ) GMOVE 
r ead ( 3 , * ) NDV , NCON 
read ( 3 , * ) IDEBUG , IPRINT 
read (3 , *) (GMAX(I) ,1=1, NCON) 

CLOSE (UNIT=3) 

OPEN (UNIT=8 , STATUS= ' OLD ' , FILE=FILNM) 

OPEN (UNIT=9 , ACCESS= ' TRANSPARENT' , FORM=' UNFORMATTED ' , FILE=FILNM1 ) 

OPEN (UNIT=1 0 , STATUS= ' OLD ' , FILE= ' HISTORY . DAT ' ) 

OPEN ( UNIT=1 1 , STATUS= ' OLD ' , FILE= ' FLAGS . DAT ' ) 

C ***************************************************************** 

C READ FILNM TO EXTRACT INFO FOR DERIVATIVES CALCULATION 
C 


READ (8 , 100) T (1) ,T(2) ,T(3) ,T(4) ,T(5) 

DO 1000 L=1 , NLINES 
C find first line of results 

IF (T ( 1) .NE. START (1) .OR.T(5) .NE. START (5) ) THEN 
READ (8, 100) T(l) ,T(2) ,T(3) ,T(4) ,T(5) 

ELSE 

READ (8, 101) 

C read some blank lines to get to beginning of data 

C initially do 10 i=l,ndv THIS SHOULD BE ACCORDING TO FILE 

if (idebug.ge. 3) PRINT *,' DESIGN VARIABLES ' 

DO 10 1=1,4 

C read the design variables X(I) 

READ (8 , 102) X(I) 
if (idebug.ge. 3) PRINT *, X(I) 

C Compute the move limits 

VLB(I) =X(I) * (1. -GMOVE/100. ) 

VUB (I) =X (I) * (1 . +GMOVE/100. ) 
if (idebug.ge. 3) PRINT *, VLB(I) ,VUB(I) 

10 CONTINUE 

C 

C ADD THE FOLLOWING LOWER BOUNDS FOR PROBLEM TO BE REALISTIC 

C 

IF (VLB (1) . LE. 2 . ) VLB ( 1 ) =2 . 0 



IF (VLB (2) .LE.0.5) VLB(2)=0.5 
IF (VLB (3) .LE.0.1) VLB(3)=0.1 
IF (VLB (4 ) . LE. 0.05) VLB(4)=0.05 

and upper limits 
IF(VUB(1) .GE.15.) VUB(1)=15.0 
IF(VUB(2) .GE.15.) VUB(2)=15.0 
IF (VUB(3) .GE. (X(2) / 2 . ) ) VUB(3)=X(2) /2 . 

IF (VUB (4) .GE. (X(l)/2.)) VUB(4)=X(1) /2 . 

Read some more blank lines 
READ (8 , 103 ) 

and then the Objective function at the design point OBJ 
and the objective at finite differences from the origin 
READ (8, 104) OBJ, (OOBJ(J) , J=1,NDV) 
if (idebug.ge. 1) THEN 

PRINT OBJECTIVE AND RESULTS OF FDs ' 

PRINT *, OBJ, (OOBJ(J) ,J=1,NDV) 

ENDIF 

convert constraints into <=0 constraints and scale 
if (idebug.ge. 1) PRINT CONSTRAINTS AND RESULTS OF FDs 

DO 11 J=1 , NCON 

READ (8, 104) G ( J) , (CONS(J,K) ,K=1,NDV) 
if (idebug.ge. 1) PRINT *, G(J) , (CONS(J,K) ,K=1,NDV) 
G(J)=G(J) /GMAX(J) -1. 

IF ( J. EQ. 2) G(J)=-G(J) 

DO 13 KK=1 , NDV 

CONS ( J, KK) =CONS ( J, KK) /GMAX ( J) -1 
IF ( J. EQ. 2 ) CONS ( J,KK)=-CONS ( J,KK) 

CONTINUE 

if (idebug.ge. 1) PRINT *, ' CORR ' ,G(J) , (CONS(J,K) 
,K=1,NDV) 

CONTINUE 

Now compute the derivatives: 

DO 12 1=1, NDV 

DF ( I ) = ( OOB J ( I ) -OBJ ) / X ( I ) / FACT 
if (idebug.ge. 1) PRINT *,'OBJ DER ',DF(I) 

DO 12 J=1 , NCON 

A ( I , J) = (CONS ( J, I) -G ( J) ) /X (I) /FACT 
if (idebug.ge. 1) PRINT DERIV ',A(I,J) 

CONTINUE 

write values to confirm 

WRITE (9) NDV, (X (I) , 1=1, NDV) , OBJ, NCON, (G(J) , J=l,NCON) , 
(DF(II) ,11=1, NDV) , ( (A(K,M) ,K=1,NDV) ,M=l,NCON) 
if (idebug.ge. 2) THEN 

print SUMMARY ' 

print *,NDV, (X(I) ,1=1, NDV) 

print *, OBJ, NCON, (G ( J) , J=l,NCON) 

print *, (DF(II) ,11=1, NDV) 

print *, ( (A(K,M) ,K=1,NDV) ,M=l,NCON) 

ENDIF 

replace values into conmin arrays and form, they will 
be passed to approx through common. 

FUNC (1) =OBJ 
DO 20 1=1, NDV 

DV ( I ) =X ( I ) 

GRAD (1,1) =DF ( I ) 



ooo 


DO 20 JJ*1,NC0N 

GRAD(I,JJ+1)*A(I,JJ) 

20 CONTINUE 

DO 21 J*l,NCON 

FUNC(J+1)=G(J) 

21 CONTINUE 
GOTO 999 

ENDIF 

1000 CONTINUE 

CLOSE (UNIT=8) 

C 

999 CONTINUE 

C INITIALIZE CONSTRAINT IDENTIFICATION VECTOR, ISC. 

DO 310 J*l,NCON+l 
310 ISC(J)«0 

C* ***************************** *************************************** 


C SOLVE OPTIMIZATION. 

350 CONTINUE 

if ( idebug . ge . 2 ) print *, 'before conmin' ,X(1) ,X(2) ,X(3) ,X(4) 

CALL CONMIN (X, VLB, VUB, G, SCAL, DF, A, S,G1,G2,B,C, ISC, IC, MSI, 

1 N1,N2,N3,N4,N5) 

if (idebug. ge. 2) print *, 'after connin' ,X(1) ,X(2) ,X(3) ,X(4) 
IF(IGOTO.EQ.O) THEN 

C reached optinum 

if ( idebug . ge . 2 ) then 

print *, 'final results' 
print *, ' ' 

print *, 'OBJECTIVE - ',OBJ 
print *,' X VECTOR ' , (X(I) ,I=1,NDV) 
print G VECTOR ' , (G(J) , J*l,NCON) 

end if 

WRITE(10, *) OBJ, (X(I) , 1=1, NDV) , (G(J) ,J=l,NCON) 
write info to new file to rerun ansys 

first, we have to read the input file for ansys and then rewrite 
it with new values 

OPEN (UNIT=4 , STATUS* 'OLD' , FILE=FILNM2) 

OPEN (UNIT=5 , STATUS* ' UNKNOWN ' , FILE=FILNM3) 

C READ AND WRITE FILE 

WRITE (5 , 110) (X(I) , 1*1, NDV) 

DO 363 NN=1 , NDV 
READ (4,*) TT 
363 CONTINUE 

DO 361 NN=1,NLINES 

READ (4 , 111 , END=362) TT 
WRITE (5, 111)TT 

361 CONTINUE 

362 CONTINUE 
ICALL=1 
REWIND (11) 

WRITE ( 1 1 , 1 1 2 ) ICALL , IFLAG 
CLOSE (UNIT=11) 



oonoo 


STOP 


ELSE 

C no convergence yet . . . 

rewind (11) 

READ ( 1 1 , 1 12 ) ICALL , IFLAG 
IF(ICALL.EQ.l) THEN 

first call to approximation, copy file and compute exponent 
IFLAG= 1 LINEAR 

2 RECIPROCAL 

3 TWO POINT EXPONENTIAL 

ICALL=0 
REWIND (11) 

WRITE ( 11 , 112) ICALL, IFLAG 
I FLAGT=I FLAG 

IF(IFLAG.EQ.3) THEN 

INQUIRE (FILE= ' SCNDGRD.DAT' , EXIST=TOF) 

IF(TOF) THEN 

OPEN (UNIT=7 , ACCESS=' TRANSPARENT' , FORM= ' UNFORMATTED ' 

1 , STATUS= ' OLD ' , FILE= ' SCNDGRD . DAT ' ) 

READ(7)NDV, (DV1(I) ,1=1, NDV) ,FUNC1(1) ,NCON, (FUNCl(J) 

1 , J=2,NCON+l) , (GRAD1 (L, 1) ,L=1,NDV) , ( (GRAD1 (K,M) 

2 , K=1,NDV) ,M=2,NCON+l) 
i f ( idebug . ge . 4 ) then 

print *, 'old point: ' , (DVl(I) ,1=1, NDV) 
print *, 'old obj. ', FUNCl(l) 
print *, 'old constr ', (FUNCl(J) ,J=2,NC0N+1) 
print *, 'old grads ' , ( (GRAD1 (K,M) ,K=1,NDV) 

1 ,M=l,NCON) 

endif 
REWIND (7) 

WRITE (7) NDV, (DV(I) ,I=1,NDV) ,FUNC(1) ,NCON, (FUNC(J) 

1 , J=2,NC0N+1) , (GRAD(L, 1) ,L=1,NDV) , ( (GRAD (K,M) 

2 , K=1,NDV) ,M=2,NC0N+1) 
if (idebug. ge. 4) then 

print *, 'Xo point: ' , (DV(I) ,1=1, NDV) 
print *, ' obj. ', FUNC(l) 
print *, ' constr ' , (FUNC(J) , J=2 ,NCON+l) 
print *, ' grads ',( (GRAD (K,M) ,K=1, NDV) 

1 ,M=1 , NCON) 

endif 

CLOSE (UNIT=7) 

ELSE 

C first call, no data in SCNDGRD.DAT yet. put it in 

IFLAG=1 

OPEN ( UNIT=7 , ACCESS= ' TRANSPARENT ' , FORM= ' UNFORMATTED ' 

1 , STATUS= ' NEW ' , FILE= ' SCNDGRD . DAT ' ) 

WRITE (7) NDV, (DV(I) ,1=1, NDV) ,FUNC(1) ,NCON, (FUNC(J) 

1 , J=2,NC0N+1) , (GRAD(L, 1) ,L=1,NDV) , ( (GRAD (K,M) 

2 , K=1,NDV) , M=2,NCON+l) 

CLOSE (UNIT=7) 

ENDIF 

ENDIF 

NFUNCS=NCON+l 
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COMPUTATION OF EXPONENT BASED ON I FLAG 

DO 710 1=1, NDV 

IF (DV(I) . EQ. 0. ) THEN 
RATX(I)=1.E8 

ELSE 

RATX ( I ) =DV1 ( I ) /DV (I) 

ENDIF 

if (idebug.ge.3) THEN 

PRINT *, 'INITIAL CALCULATIONS IFLAG= ' , I FLAG, 'X(I) = ' 
,X(I), I, 'DV1 (I) /DV(I) ' ,RATX(I) 

ENDIF 

IF( (IFLAG.EQ.l) .OR. (X(I) .EQ.O.) .OR. (RATX(I) .EQ.l.) ) 

THEN 

if (idebug.ge.2) print In linear code ' 

LINEAR APPROXIMATION 
DO 711 J=1 , NFUNCS 
P(I,J)=1. 

CONTINUE 

ELSE 

IF( (IFLAG.EQ.2) .OR. (DV(I) .EQ.O.)) THEN 

if (idebug.ge.2 ) print *,' In Reciprocal code ' 
RECIPROCAL APPROXIMATION 
DO 712 J=l, NFUNCS 
P(I,J)=-1. 

CONTINUE 

ELSE 

if (idebug.ge.2 ) print In 2 point code ' 

2 POINT EXPONENTIAL APPROXIMATION 
DO 713 J=l, NFUNCS 

IF(GRAD(I, J) .EQ.O.) THEN 
P(I,J)-1. 

ELSE 

RATDER ( I , J) =GRAD1 ( I , J) /GRAD ( I , J) 

IF ( (RATX (I) .LE.O.) .OR. ( RATDER ( I, J) .LE.O.) ) 
THEN 

P(I # J)-1. 

ELSE 

P(I , J) =DLOG (RATDER (I ,J) ) /DLOG (RATX (I) ) +1 
IF(P(I, J) .GE.l.) THEN 
P(I,J)=1. 

ELSE 

IF(P(I,J) .LE.-l.) P(I,J)=-1. 

ENDIF 

ENDIF 

ENDIF 

CONTINUE 

ENDIF 

ENDIF 

if (idebug.ge.2 ) PRINT *,'I, EXPONENT ****** ,I,P(I, l) 
CONTINUE 


IFLAG=IFLAGT 

ENDIF 

IF (INFO.EQ. 1) THEN 
AF (1) =FUNC (1) 

if (idebug.ge.3) PRINT*, 'OBJ ',obj 
DO 359 J=l,NCON 



AF ( J+l) =FUNC ( J+l) 

if ( idebug . ge . 3 ) PRINT* , ' CONS #',j, G(J) 

359 CONTINUE 

if (idebug. ge. 2) PRINT *,'CALL TO APPROXIMATION ' 

C this is the call to the approximation 

CALL APPROX (X,AF,P, NDV, NCON) 

C Resubstituting values in OBJ and CONS 

OBJ=AF(l) 

if (idebug. ge. 3) PRINT *,'OBJ ',obj 
DO 360 J=l,NCON 
G(J)=AF(J+1) 

if (idebug. ge. 3) PRINT *,'CONS #',j, G(J) 

360 CONTINUE 
ELSE 

if (idebug. ge. 3) PRINT *, '# Info ne 1 ??? ',INFO 
ENDIF 
ENDIF 
GOTO 350 

C FORMATS 

98 FORMAT (F3.0) 

99 FORMAT (A12/A12/A12/A12) 

100 FORMAT (5A4) 

101 FORMAT (IX,//) 

102 FORMAT (5X,E12. 6) 

103 FORMAT (IX,/////////) 

104 FORMAT (5X,6E13. 6) 

110 FORMAT ('H=' ,E12.6/'B=' , E12 . 6/ 'T=' , E12 . 6/ 'D=' , E12 . 6) 

C110 FORMAT('H=',E12.6/'B=',E12.6) 

111 FORMAT (A80) 

112 FORMAT (2 12) 

C 

END 

SUBROUTINE APPROX ( AV , AF , P , NDV , NCON) 

C 

C THIS SUBROUTINE IS CALLED FROM OPTRUN TO PERFORM 

C VARIOUS APPROXIMATIONS OF THE FUNCTIONS (OBJECTIVES 

C AND CONSTRAINTS) . A FLAG WILL SELECT LINEAR, 

C RECIPROCAL OR IMPROVED APPROXIMATION. TWO SETS OF DATA 

C ARE NEEDED SINCE THE IMPROVED APPROXIMATION RELIES ON 

C PAST ANALYSES TO IMPROVE THE APPROXIMATION. 

C 

C Georges Fadel June 1989 

C Oct 1990 

C Jan 1991 

C AV is the vector of VARIABLES 
C 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DIMENSION AV (4) , AF(7), P(4,7) 

COMMON /INFOIN/ DV(4) , FUNC(7) , GRAD(4,7) 

COMMON /FLAGS/ IFLAG, ICALL, IDEBUG 
NFUNCS=NCON+l 

C 

C NOW THE EXPONENT IS KNOWN, LETS COMPUTE THE APPROXIMATING FUNCTION 
C 

if (idebug. ge.l) PRINT*, ' EXPONENT KNOWN ' 

DO 30 J=1 , NFUNCS 
DO 20 1=1, NDV 

IF ( (P (I , J) . EQ . 1 . ) .OR. (ABS (P (I , J) ) .LE. 0. 00001) ) THEN 
AF ( J) =AF ( J) + (AV(I) -DV(I) ) *GRAD(I, J) 

ELSE 



IF(P(I, J) .EQ.-l.) THEN 

AF(J)=AF(J)+(AV(I)-DV(I))*(DV(I)/AV(I) )*GRAD(I, J) 
ELSE 

AF(J)=AF(J)+((AV(I)/DV(I))**P(I,J)-1.) 

1 *DV(I)*GRAD(I,J)/P(I,J) 

ENDIF 

ENDIF 

20 CONTINUE 
30 CONTINUE 
999 RETURN 
END 



APPENDIX C 
OPTIM.DAT file 


BEAM4 . OUT 
BEAM 4 . INP 
BEAM4 . OLD 
BEAM4 . DAT 
50 . 

4 

3 

0 4 

5.0 30 . 80 . 



APPENDIX D 

OPTIMI.BAT batch file to execute optimization 


echo off 
els 

echo OPTIMIZATION BATCH FILE TO TEST APPROXIMATIONS WITH ANSYS AND CONMI 
echo £ 


echo START OF OPTIMIZATION PROCEDURE 

echo £ 

echo call to ANSYS with initial design variables set in file 

echo xxxxxxx.dat in ansys format 


ERASE %1.0UT 
echo ^ 

call ANSYS -I %l.dat -O %l.out 
COPY %1.DAT % 1 . OLD 
echo £ 

echo First ANSYS run COMPLETED. Results are written to %l.out 

echo £ 

: loopl 

call browse %l.out 
echo £ 

echo +++++++++++++ IN LOOP +++++++++++++++++ 

echo £ 

echo call optimization program using design variables and derivatives 

echo £ 

call optrun2 

echo £ 

echo £ 

echo — — — CONVERGED IN APPROXIMATION LOOP 

copy hist.dat+history.dat hist.dat 
echo £ 

ERASE %1.0UT 

call ANSYS -I %1.DAT -0 %l.out 
echo £ 
echo 
echo £ 

REM if not converged 
goto loopl 
REM else 
stop 


RERUN ANALYSIS (ANSYS) 



